INSTITUT NATIONAL DE RECHERCHE EN INFORMATIQUE ET EN AUTOMATIQUE 

The Multidimensional Refinement Indicators 
Algorithm for Optimal Parameterization 

Hend Ben Ameur — Francois Clement — Pierre Weis — Guy Chavent 



N° 5940 

Juin 2006 
Themes NUM et SYM 



0X# 1 J\ Kiri 

^H^B^ ROCOUENCOURT 



The Multidimensional Refinement Indicators Algorithm for Optimal 

Parameterization 

Hend Ben AnieuiB , Francois Clemen!! , Pierre Weis|,Guy ChavenH 

Themes NUM et SYM — Systemes numeriques et Systemes symboliques 
Projets Estime et Cristal 

Rapport de recherche n° 5940 — Juin 2006 — [TBlpages 



Abstract: 

The estimation of distributed parameters in a partial differential equation (PDE) from measures of the solution of the 
PDE may lead to underdetermination problems. The choice of a parameterization is a frequently used way of adding 
a priori information by reducing the number of unknowns according to the physics of the problem. The refinement 
indicators algorithm provides a fruitful adaptive parameterization technique that parsimoniously opens the degrees of 
freedom in an iterative way. We present a new general form of the refinement indicators algorithm that is applicable to the 
estimation of distributed multidimensional parameters in any PDE. In the linear case, we state the relationship between 
the refinement indicator and the decrease of the usual least-squares data misfit objective function. We give numerical 
results in the simple case of the identity model, and this application reveals the refinement indicators algorithm as an 
image segmentation technique. 
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L'algorithme des indicateurs de raffinement multidimensionel pour une 

parametrisation optimale 

Resume : 

L' estimation de parametres distribues dans des equations aux derivees partielles (EDP) a partir de mesures de la 
solution de FEDP peut mener a des problemes de sous-determination. Le choix d'une parametrisation est un moyen usuel 
pour ajouter de F information a priori en reduisant le nombre d'inconnues en relation avec la physique du probleme. 
L'algorithme des indicateurs de raffinement fourni une technique de parametrisation adaptative fructueuse qui ouvre 
parcimonieusement les degres de liberte de facon iterative. Nous presentons une nouvelle forme generate de l'algorithme 
des indicateurs de raffinement qui s'applique a l'estimation des parametres multi-dimensionnels dans toute EDP. Dans le 
cas lineaire, nous etablissons le lien entre l'indicateur de raffinement et la decroissance de la fonction objectif des moindres 
carres quantifiant Ferreur aux donnees. Nous donnons des resultats numeriques pour le cas simple du modele identite, 
et cette application permet de voir V algorithme des indicateurs de raffinement comme une technique de segmentation 
d' image. 

Mots-cles : probleme inverse, parametrisation, raffinement optimal, segmentation d' image, programmation fonctionnelle 
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1 Introduction 

The inverse problems of parameter estimation consist in identifying unknown parameters which are space dependent 
coefficients in a system of partial differential equations modeling a physical problem. The parameter estimation problem 
can be set as a minimization problem of an objective function defined as the least-squares misfit between measurements 
and the corresponding quantities computed with a chosen parameterization of the parameters (solution of the partial 
differential equations system under consideration). 

Experimental measurements are expensive, thus one of the difficulties when solving a parameter estimation inverse 
problem is that the data is usually insufficient to estimate the value of the parameter in each cell of the computational 
mesh, the parameter estimation problem is underdetermined. Therefore, a parameterization of the sought parameter is 
chosen to reduce the number of unknowns. The multiscale parameterization approach is popular for many applications, 
e.g. see lfT5l l6l [TTl [D. It consists in solving the problem through successive approximations by refining uniformly the 
parameter at the next finer scale all over the domain and stopping the process when the refinement does not induce 
significant decrease of the objective function any more. This method is very attractive as it may provide a regularizing 
effect on the resolution of the inverse problem and it may also circumvent local minima problems as shown in ifTolfTTI . 
But its main drawback is the over-parameterization problem due to the uniform refinement of the parameter, see Q [8) . 

To avoid this pitfall, the refinement indicators algorithm provides an adaptative parameterization technique that parsi- 
moniously opens the degrees of freedom in an iterative way driven at first order by the model to locate the discontinuities 
of the sought parameter. The detailed definition of the technique was described in |2 ] for the application to the estimation 
of hydraulic transmissivities, and a variant of the idea was first briefly presented in J9l for the application to the history 
matching of an oil reservoir. 

In the present work, we extend the definition of refinement indicators to the more general case of multidimensional 
distributed parameters, we make the link between the (first order) refinement indicator and the exact decrease of the 
objective function in the linear case, and we propose a generic version of the algorithm that may be applied to any inverse 
problem of parameter estimation in a system of partial differential equations. We show numerical results for the simple 
case of the identity model. The application of the refinement indicators algorithm to the identity model applied to images 
is interpreted as an amazing image segmentation technique. 

The paper is organized as follows. Sectionals devoted to theoretical developments about adaptive parameterization 
through the refinement indicators notion: we quantify the effect of refinement on the decrease of the objective function, 
and we generalize the definition of refinement indicators for multidimensional parameters. We present in Section [3] the 
general form of the refinement indicators algorithm, we discuss the issue of finding the best cutting for multidimensional 
parameters and briefly point out implementation subtleties of the algorithm. In section|4] we focus on the identity model; 
we then apply the multidimensional refinement indicators algorithm to the particular case of RGB color images. 

2 Adaptive inversion and refinement indicators 

The inverse problem of parameter estimation is first set as the minimization of a least-squares objective function. Then, 
the notion of adaptive parameterization induced by zonations of the domain is presented in the discrete case for piecewise 
constant parameters. The decrease of the optimal value of the objective function resulting from the splitting of one zone 
into two subzones is then quantified for the linearized problem in the continuous case, and refinement indicators are finally 
defined in the general case from a first order approximation. 

2.1 Least-squares inversion 

When defining a system of partial differential equations (PDE) modeling a physical behavior in a domain fict" (e.g. 
with N = 1, 2 or 3), one usually introduces parameters such as material properties. We consider the case where these 
parameters are distributed and possibly vector valued, i.e. belong to a space P of functions defined over £2 with values in 
W'p (n p > 1 is the dimension of the vector parameter p(x) at any x e £2). 

Some of these parameters may not be accessible by direct measurements, so they have to be estimated indirectly. 

If 

Ptrue £ P denotes such an unknown parameter that we are looking for in a set of admissible parameters P ad , and if 
d ~ J (/?true) 6 (the data) denotes the corresponding vector of available measurements of the solution of the PDE, 
one can attempt an indirect determination of p tme from d by solving the least-squares inverse problem: 

(1) minimize J(p) for p € P ad 
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where J(p) is the least-squares misfit between d and the corresponding quantities J (p) computed from the current pa- 
rameter p, 

(2) 3(p) = \\\d-T(p)\\ 2 - 

The forward operator J is the composition of the model operator M (which computes the solution of the PDE for 
a given parameter p), with the observation operator O (which computes the output of the observation device applied to 
the solution of the PDE). For example, observations can be made of the values of the solution of the PDE at a set of 
measurement points of the domain £2. 

Minimizing such an objective function avoids the difficulty of inverting £W , which is impossible in most cases. 

2.2 Adaptive parameterization 

In parameter estimation problems, the data is usually insufficient to estimate the value of the (possibly vector valued) 
unknown parameter at each point x G £2 (i.e. in each grid cell after discretization). It is usually impossible to increase the 
number of measurements, and one solution is then to reduce the number of unknowns by searching for the parameter in a 
subspace of P of finite (and small) dimension. 

For that, we proceed in two steps: first we construct a finite dimensional subspace P„ of the infinite dimensional 
space P, then we look for an approximation p„ of the unknown and infinite dimensional true parameter in P„ d — P„ nf ad . 
The dimension of P n should remain small enough to avoid over-parameterization, i.e. the Jacobian of J as a function 
from P n to M" d should have a full rank, and the approximation p„ should explain the data up to the noise level. Notice that 
here both P„ and p„ are unknown. 

For a given subspace P„, a natural choice for p n is a minimizer of the objective function J in P* d . 

It is classical to build the subspace P n in an iterative way by successive refinements as in the multiscale approach, 
see [15 6| El Q] • We use the adaptive parameterization technique developed in J2] . These approaches are also known to 
have a regularizing effect on the resolution of the inverse problem, e.g. see |[Kfl[l31[H^[T4"l l5l. 

The adaptive parameterization approach consists in adding one (vectorial) degree of freedom at a time. The new 
degree of freedom is chosen according to refinement indicators which evaluate the benefit for the objective function of 
the introduction of a given family of degrees of freedom. Hence, each subspace P„ is of dimension n x n p . The process is 
stopped when p n explains the data up to the noise level. 

It is convenient to consider P„ as the range of an — unknown — parameterization map 

(3) $ n : m„ G Mf .— > p n G Pf 

where m„ is the coarse parameter (of small finite dimension), by opposition to the fine parameter p„ (of large, and possibly 
infinite dimension). M^ d is the space of admissible coarse parameters. Typically, m„ is made of the coefficients of p n G P n 
on a basis of P n , in which case fP„ is a linear operator. But the parameterization map can also be nonlinear, as for 
example in the case where m„ is made of the coefficients of a closed form formula defining the fine parameter p„. For any 
parameterization map !P n , we define the same objective function on Mff by 

(4) 3 n {m n )=3{fE n {m n )), 
and the least-squares problem becomes: 

(5) minimize J„(m n ) for m„ G Mj} d . 

Let Z„ = (Z„j)i<j<„ be a partition of the closure of the domain £2 in n parts, i.e. for all j G {1, . . . ,«} Z n j C £2, 
\J l< j <n Z n j = £2, and Z„ j fiZ^ = as soon as j ^ k. We consider that the subspace P n is made of piecewise constant 
(vector valued) functions on this partition. We call zonation the partition Z„, and zone each part Z„ j of the zonation. The 
parameterization map associated with the zonation Z„ is then 

(6) <B n : m n — (m n j)i<j<„ i — ► p„ — (x G £2 i-> m n j G M" p when x G Z„j). 

It associates the coarse parameter m„ with the function p n which takes on each zone Z nj the constant value m„j. 

In the iterative process, when going from P„ to P n +i, we introduce a new (vectorial) degree of freedom by dividing 
one zone of Z„ into two subzones, thus producing a new zonation Z n+ \ having one zone more than Z n . 
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2.3 Quantifying the effect of refinement for the linearized problem 

Considering a current zonation, the best refinement would be the one corresponding to the largest decrease of the objective 
function. This decrease can only be computed by actually solving the least-squares problem (0 for each possible refine- 
ment to be tested, hence only a very small number of them could be evaluated before reaching a decision. So, we search 
for a closed form formula for the decrease of the objective function for the linearized problem. We shall still denote by J 
the linearized forward map. 

For simplicity, we consider the case where the optimization problem © is unconstrained, so that P* d = P n and Mff = 

R nxn p for a jj n 

Let us compute the decrease of the optimal objective function (0]i resulting from refining the current zonation by 
splitting one of the zones into two subzones by means of a continuum of intermediate minimization problems constraining 
the discontinuity jump at the interface between the subzones. 




m\ on Z\ 





mn on Z22 



TO21 onZ2i 




Figure 1 : Refinement of a single zone zonation represented in a two-dimensional space for simplicity. Left: one-zone 
parameterization <2\ with Q, = Z\\. Right: two-zone parameterization ¥2 with Q, = Zi^\ UZ2,2- 



Let us consider the case where only one zone covers the whole domain, corresponding to the one-zone parameteriza- 
tion <2\, see Figure |231 (left). The refinement builds the two-zone parameterization 2>2 of Figure |231 (right). We denote 



by m° pt : 



(in 



opt 



1,1) md h 



opt opt 



(m^lim^) 7 an d J^ 1 " tne optimal coarse parameter and objective function respectively 



when considering the parameterizations i>i and 2>2. 



If the discontinuity jump c opt = m 



opt 

2,1 



opt 
' 7 2.2 



was known, then minimizing J2 (i.e. minimizing J considering the 



parameterization T2) under the constraint n%2,\ — /«2,2 = c opt would give us the same optimal coarse parameter m? p and the 
same optimal objective function value J^, the ones obtained without any constraints. And minimizing J2 and under the 
constraint W2. 1 — W2,2 = would keep the optimal values m ° pt and J7° pt , obtained with the parameterization ¥\ . Thus, when 
the discontinuity jump c goes continuously from to c opt , then the minimum of J2 under the constraint W2. 1 — /W2.2 — c 8 oes 
continuously (and even in a continuously differentiable way when J is continuously differentiable) from the minimum 
obtained for a single zone to the minimum obtained with the two zones. 

The components 1112.1 and m.2,2 of the coarse parameter are column vectors of dimension n p (i.e. n p = 1 in the scalar 
case). We can denote the constraint on the discontinuity in matrix form by 



(7) 



Am2 



where m2 — {m-2,1, m2,2) T and the n p x 2n p rectangular matrix A is of the form 



(8) 



A = 



( 1 


Vo 



-1 
: 



\ 



■■ 

1 









-1 / 



We denote the (linear) direct operator in terms of the coarse parameter m„ by J n = J o <2 n . In the following computa- 
tion of the decrease of the optimal value of the objective function, we will only consider the two-zone parameterization i>2- 
Hence, when no ambiguity can arise, we will simply write f for J2, 3 for 32 and m for TO2. The gradient of the quadratic 
objective function J (— J2) is given by 

(9) Vy(m) = J T {jm-d). 

The Lagrangian function associated with the minimization of J7 under the constraint (Q writes 

(10) L c (m,X) = j(m) + (X,Am-c) 
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where A, is the Lagrange multiplier associated with the constraint. It is a vector of dimension n p . Then, the Lagrange 
condition ensures that at the optimum, m c opt = argminAm=c 3 (fit) and A.^ opt are obtained by solving 



(11) 



— (m,X) = Vj/ (m) +A T X = 0, 
dm 

dL c 

-zrr- {m, X) = Am — c = 0. 
dX 



Using the expression (0 of the gradient, the Lagrange condition for c = becomes 



(12) 



J T (T m°- opt - d) +A T X°' opt = 0, 
Am !Op t = Q 



The last equation writes m^ 1 — iru' < ? t (= m i P i)> and we have 

(13) 3° opt = 3. F2 (m Qopt ) = 3l 



The optimal coarse parameter m opt (= ni^ 1 ) = argmin J (m) (i.e., without constraint) satisfies the optimality condition 

(14) Vj7(m opt ) = f T (Tm o ' 1,t -d)=0. 
Let e = m opt — m opt . Developing the squared norm yields 

3 opt = 3 ' opt +^\\fe\\ 2 + (fm > o v t -d,!Fe). 

But, from the optimality condition (fT4l) . we have 

(fm°-°P t -d,Te) = -\\!re\\ 2 
and taking the difference between ( TT~4T > and the first equation of dT2b leads to 

Therefore, when f T fis invertibl^lL we finally obtain 

(15) 3T-3T=\\\^™%T? r , = \{A T X^ 

where the norm and the scalar product are now defined in the space of coarse parameters for the refined zonation (and not 
in the space of data as before). 

Remark 1 In the general case, an n-zone parameterization is refined into an (n + \ )-zone parameterization by splitting 
only one zone into two subzones. Then, the matrix defining the new discontinuity jump of the parameter is a block matrix 
containing only one nonzero block corresponding to the discontinuity of the parameter between the two subzones and this 
block is equal to the matrix A o/dS]). Hence, the computation of the decrease of the optimal value of the objective function 
is exactly the same. 



Because the number of zones is supposed to stay small, after discretization, the direct operator J n is usually injective. 
But applying (ij Tn) 1 is at least as expensive as solving the linearized minimization problem. Thus, the exact compu- 
tation of the decrease of the optimal value of the objective function for the linearized problem is generally not compatible 
with the idea of a very fast computation of an indicator on the quality of the new zonation. Nevertheless, an important 
exception is the case of the identity direct operator; it will be developed in section|U The next subsection is devoted to a 
first order computation in the general nonlinear case. 

'E.g., when f is injective in the finite-dimensional case. 
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2.4 Refinement indicators for multidimensional parameters 

Following the notations of the previous subsection, let J e opt denote the optimal value of the objective function obtained 
with the parameterization r?2 under the constraint (0. A first order development of j c o p 1 with respect to the discontinuity 
jump c at c = can be written 

3 7 e,opt 

J copt = _ 7 0,opt + ^ — c + ><< 
oc |c=0 

The norm of the quantity — ^ — tells us, at the first order, how large would be the difference between y c '°P l 

dc \ c =o 

and J7 0,opt = .7° pt (from (fTSTl). So, it gives us the first order effect on the optimal value of the objective function produced 
by refining and allowing a discontinuity of the parameter between the two subzones of the parameterization Similarly 
to the particular scalar case developed in [2], this norm is called the refinement indicator associated with the splitting of 
the zone Zu into the two subzones Z2.1 and Z2.2, and it is denoted by / . 
Deriving the expression (TTOb of the Lagrangian with respect to c gives 

d£C ( \\ dj l \ \ 
— — (to, A,) — — (to) — k. 

oc oc 

Hence, from the Lagrange condition (fTTT i. we have at the optimum for c = 



Then, since 



- -(m^-X^^O. 
oc 



dc | c =0 dc \c=0 dc 

this shows that the refinement indicator is nothing but the norm of the Lagrange multiplier, 

(16) I = \\k°<*\\. 

Remark 2 In this section, we have considered that the vector parameter was a piecewise constant function (i.e. constant 
in each zone). The idea of refinement indicators is based on defining a zonation by identifying the interfaces between the 
zones that are related to the discontinuities of the parameter. To estimate more regular parameters, one can use similar 
refinement indicators in the more general case where the vector parameter is a piecewise higher order function, which is 
continuous in each zone and presents discontinuities at the interfaces between the zones. 

3 Refinement indicators algorithm 

To solve numerically the parameter estimation problem, we consider a discretization of the domain £2 by a fine mesh 
Tf, = (Ki)i£i, i.e. with £2 = U;e/^!- Let nj = card/ be the number of cells. The fine parameter is then the vector p = (pi)i£i 
which approximates the parameter on each cell Kj of T/,. We consider zonations Z„ = (Z„j) i</<« following the mesh such 
that the associated parameterization map (0) is 

<B n : (R"") n — > (Wp)"' 

(17) 

m„ = (m n j)i<j<„ 1 — > p n = (p n ,i)iei with p„j = m n j when Kj C Z„j. 
In other words, when dealing with a discrete problem, a zonation is related to a partition of the set / of indices: 

Z„j = [J K{ with I = LJ Ifi,j (disjoint union). 

iel nJ l<j<n 

The finite set of measurement points is indexed by /,„ C /. Keeping the same notations, the discrete objective function 
to minimize writes 

(is) Hp) = ^ 1 £(di-ui) 2 

L ier m 

where di is the measurement in the cell Kj and Uj denotes the discrete approximation of J (p) in the same cell Kj. 
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3.1 The adaptive parameterization technique 

In practice, the optimal coarse parameter ra° pt associated with the current parameterization fZ>„ is computed by applying 
a gradient algorithm to the minimization of the corresponding objective function J7 M . Hence, not only m° pt = ('w° P j)i<;<n 
is available, but also the coarse gradient Vj7„(ra° pt ), which vanishes at the reached minimum. The key is to compute the 
gradient of J n for the current parameterization by the adjoint approach considering the fine discretization 1), of £2. This 
approach provides, at no additional cost, the fine gradient 

(19) Vpj[pT)= (j^, {pT) ). eI 

where p° pt = £P„m° pt , since the coarse gradient of J7„ is given by Vj„(m„) = V^Vj (Pn), and hence, it is simply obtained 
by summing the components of the fine gradient inside each zone of the current parameterization £P„. 

A refinement of the current parameterization <E n is obtained by cutting one zone Z„j into two subzones Z™j and Z™j 9 
with 1 < /'i , }% < n + 1. The resulting (« + 1 )-zone zonation Z™' is a candidate for Z„ + i . This operation is fully character- 
ized by the subsets of indices corresponding to the new subzones which satisfy the disjoint union I„j = iffi Ll/Jfjj. The 
associated parameterization is denoted by V™ 1 and the corresponding objective function by J7 n cut . Opening the new (vec- 
torial) degree of freedom while keeping the same value in the two subzones does not change the optimum of the objective 
function, and the components of the gradient Vj/ n cut ((m™ t ) opt ) are all equal to zero except for components j\ and ji that 
are opposite vectors (of dimension n p ). From the first equation in ( fTTT i and the form of matrix A, we deduce that the 
Lagrange multiplier (A.™ 1 ) 0,01 " is equal to the component 72 of the previous gradient. Hence, the refinement indicator J° ut 
associated with the current refinement can be computed easily, at almost no additional cost, from the partial derivatives of 
the objective function with respect to the fine parameter in the subzones by 

(20) C' = 

Thus, we can define a large number of tentative cuttings producing possible refinements in each zone, and compute their 
corresponding refinement indicators. 

Since the refinement indicators give us only a first order information on the decrease of the optimal objective function, 
one can select a set of cuttings associated with some of the highest values of the indicators, and not only with the largest 
one. Then, the objective function is minimized for all the parameterizations obtained by implementing each of these 
selected cuttings. And finally, the cutting defining the next zonation is the one leading to the largest decrease of the 
objective function. In the general case, this minimization phase is very expensive. It is thus important to keep the number 
of selected cuttings very small. 

The initial value of the coarse parameter for these minimizations is obtained from the previous optimal value by 
duplicating the value in the split zone, i.e. 

/ cut \init / cut \init opt „„ j / cut \init opt f ./ / . 

(21) \ m n,h) = ( m n,j 2 ) = m nj and ( W „J'J = m »J for J T Jl ,J2- 

The parameterization given by ( fTTb is injective, and we can define its least-squares pseudoinverse by T n = !P^)~ l . 
It computes the mean value on each zone. When considering the scalar product weighted by the measure of each cell of 
the fine mesh, the previous pseudoinverse corresponds to the projection map given by 



(p7) 



dpi 



ipT) 



i n : (R"p )"' — ► (R"p) n 

(22) , N , ., Liei n] meas i K i)Pn,i 

Pn = [Pn,i)iei 1 — > m„ = {m n j)i<j<„ with m n j = '■ — — . 

meas(Z n; j 

Then, OTT) rewrites in the more compact form (mj; ut ) lnlt = fp^ ut p^ pt . 

At the end, the next (n + l)-zone parameterization map fp„ + i has only one (vectorial) degree of freedom more, but it 
produces a significant benefit effect on the objective function. 



Remark 3 The adaptive parameterization technique is independent of the fine discretization T/, of the domain £2. In 
particular, it is valid for regular or irregular meshes, structured or unstructured meshes. 
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3.2 A generic algorithm for adaptive parameterization 

The following algorithm can be applied to the estimation of distributed multidimensional parameters in any partial differ- 
ential equation. 

Initialization 

-1 . Choose an initial no-zone parameterization iP„ and an initial coarse parameter m™'. 

0. Compute ra° pt = argmin j„ (m„ ) from m™ 1 and pjjj* = 2n m°f. 
Iterations For n > no, do until convergence: 

1. Compute the fine gradient g° p = V P J (p° pt ). 

2. Choose a set of cuttings C„ defining parameterization candidates (&„)keC n - 

3. For all k 6 C„, compute the corresponding indicator J* from g° pt using ( f20b . 

4. Compute the largest indicator / n max . 

5. Select a subset C„(/„ max ) C C„ of cuttings associated with the highest indicators. 

6. ForallfeeC„(/„ max ), 

compute {m k n )° ?x — argmin J*(m*) from (m*) lnlt = i k p°? 1 - 

7. Keep the best cutting ko = ^E m ^ n keC„(i ms!L ) ■?n(( m n)° pt ) 
and set fP n+1 = m*, = (m* ) *" and = fP^m^. 

This generic algorithm may be adapted by the user to match his/her needs. The convergence criterion can be specified, 
for instance: a small value of the objective function (meaning that the data are correctly fitted), a small value of the largest 
refinement indicator (meaning that there is no more important discontinuities left to discover), the maximum number of 
zones is reached. . . In addition, the user is free to choose the refinement strategy defining the full set of cuttings and the 
subset of those with the highest indicators. 

Furthermore, as in [2], it is possible to add at the end of each iteration a coarsening step that allows zones corresponding 
to coarse parameters with close values to be merged. 

3.3 Best cutting for multidimensional parameters 

Let us suppose that we want to find among all possible cuttings the one corresponding to the best indicator. For example, 
this is useful in section [4] with the identity direct operator for which this best cutting is related to the best decrease of the 
objective function. But this absolute best cutting may also be useful in the general case to quantify the adequacy of the 
chosen set of cuttings to the inverse problem trough the ratio of the best indicator within the chosen set of cuttings to the 
absolute best indicator. A high ratio shows that there exists cuttings outside the chosen set that provide higher indicators. 

In the scalar case, it is obvious from equation ( f20b that the best cutting for a given zone corresponds to follow the sign 
of the fine gradient (fT9] i. In the multidimensional case, this is no more so simple as finding the best cutting amounts to 
solve a discrete optimization problem. We give here several possible heuristics. 

The simpler idea consists in defining the sign of a vector and then to use the same technique as in the scalar case. 
There are several possible definitions for the sign of a vector. Of course, none of them are fully satisfactory, but they can 
provide pertinent zonations, as it can be seen in section [4] where the sign of a three-dimensional vector is obtained by 
majority vote of the signs of all components. It could also have been defined as the sign of the sum of all components of 
the vector. 

Another simple idea is to consider a zonation for each component of the parameter. Then, the algorithm for an 
n^-dimensional parameter on a grid of size nj would be equivalent to the algorithm for a scalar parameter on a grid of 
size n p x «/. Again, with the example of section 2] where the data are colored images, n p = 3 and nj is the number of 
pixels. In this case, the algorithm for an RGB (Red/Green/Blue) image of size nj would be equivalent to the algorithm for 
3 grayscale images of the same size n/. This kind of decomposition reminds the dawn of color photography where three 
black-and-white pictures were taken with three colored filters, and more recently the 3-LCD video projectors. The main 
drawback of this technique is that the compound zonation corresponding to the superimposition of the zonations for all 
the components may contain a lot of small zones. Indeed, the intersection of n p zonations with n zones each can produce 

up to n"p zones. In other words, we would add up to Ylk'=o C k p n k degrees of freedom per component at each iteration! 
Which is completely in contradiction with the adaptive parameterization goal of providing the coarsest parameter, i.e. 
adding at most only one degree of freedom per component per iteration. 
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Keeping the attractive idea of computing a local indicator for each — scalar — component of the parameter separately 
we can select at each iteration among the best cuttings for each component the one providing the highest global indicator. 
Then, this cutting is applied to all components, thus guaranteeing the addition of a single degree of freedom per compo- 
nent. Of course, this heuristics does not always provide the highest possible global indicator, but it seems reasonable to 
us and we suggest to use it. 



3.4 Implementation issues 

The first issue is the genericity of the implementation of the algorithm. Numerical results shown in [2 1 were obtained with 
a Fortran 77 program where the refinement algorithm was hard-coded within the model for the very specific problem of 
the estimation of hydraulic transmissivities in a parabolic flow equation set on a planar regular mesh. Once the technique 
has been validated on this example, it is very advisable to propose to the community a generic implementation that will 
allow to apply the refinement indicator adaptive parameterization algorithm to any problem implemented in a wide variety 
of programming languages. 

The second issue is related to the definition of the set of cuttings, and more particularly to the implementation of 
it. In the general case, it is not possible to compute the true parameterization, and then, the adaptive parameterization 
technique is based on the definition of a set of cuttings, which would preferably be adapted to the problem. In (2), 
the two-dimensional mesh was supposed regular and rectangular, and the refinements were of four elementary kinds: 
horizontal cutting, vertical cutting, oblique cutting and checkerboard cutting. At each iteration, the indicators related to 
all possible cuttings within the selected family for all current zones were computed. Testing larger sets of cuttings, e.g. by 
enriching the elementary cutting offer, is very interesting because it will allow the user to use a priori information adapted 
to the problem. Furthermore, the ability to specify that some elementary cuttings should only apply to a particular part of 
the domain will allow to interpret more closely the a priori information. 

The last issue concerns the performance of the implemented code. At each iteration, steps [3] and [6] of the algorithm 
correspond respectively to the computation of all the indicators associated with the chosen cuttings and to the actual 
minimization for the selected cuttings. Both steps involve fully independent computations that may run in parallel. 

All these aspects show the necessity for high-level programming capabilities. In [3], we choose Objective Caml 
(OCaml) which is a fast modern high-level functional programming language based on sound theoretical foundation 
and formal semantic^ It is particularly well suited for the implementation of complex algorithms. The generic driver is 
written in OCaml and data exchange with the worker follow a simple and safe communication protocol that is imple- 
mented in several common programming languages including OCAML, C, C++ and FORTRAN. Here the generic tasks 
sent to the worker are "cost", "grad" and "optim" respectively to implement the functions p i — ► J (p), p i — > Vj7 (p) 
and (m,fp) i — > (m opt ,j opt ) = (argminj7j> (m), J? (m opt )). The design of a mini-language allow a flexible and convenient 
way to define the set of cuttings. Furthermore, automatic parallelization capabilities through the skeleton programming 
system provided by OCamlP3l0 has recently been successfully used in the field of Scientific Computation for domain 
decomposition applications, see lfT2l[T3l . 



4 Application to the identity model 

To illustrate the use of our OCAML implementation, detailed in [3|, and the capabilities of the refinement indicators 
algorithm in the multidimensional case, presented in the previous section, we consider the simplest model ever: the 
identity model. Thus, in this section, we assume that J = Id. Of course, in this very simple case, the inverse problem does 
not present any difficulty as the minimum of the objective function (fl8T l is exactly zero, and it is obtained for p° p{ = d. Note 
that the fine gradient is then simply given by the difference Vj7 (p) =p — d. Nevertheless, assuming that the data is actually 
an image of size nj = n x x n y , then the mesh is regular and rectangular, the cells are called pixels, and the application of 
the refinement indicators algorithm amounts in this case to determine groups of pixels with the same color — or gray level. 
This is the basic idea of image segmentation for edge detection, shape recognition, compression/simplification, and so 
on. . . 

The segmentation of a grayscale image is simply a scalar example, but even in this case, the determination of the 
best cutting providing the highest decrease of the objective function is not an easy task, see J4J. The case of color 
images is more complicated as the color information is typically multidimensional. In computer applications, it is usually 
represented using a three-dimensional color space — hence n p = 3: either RGB, for the Red/Green/Blue color cube, HSV, 
for the Hue/Saturation/Value color cone, or HSL, for the Hue/Saturation/Lightness color double cone. Of course, it is 

2 For which the best cutting associated with the highest local indicator simply follows the sign of the — scalar — fine gradient for each component. 

3 See http: / /www.ocaml .org/. 
4 See http : / /www . ocamlp31 . org/ . 
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possible to define a scalar discrete representation of colors, e.g. by packing the 24-bit RGB encoding into a single integer 
ranging from to 2 24 — 1 . But, such a nonlinear mapping produces weird effects when averaging on a zone: the mean 
color between red pixels and blue pixels may become an odd green, instead of the natural purple value. We consider here 
the RGB representation of color images. 

In the case of the identity model, the operator fjj n of subsection 12.31 reduces to the diagonal matrix (P„^P n which 
gives the numbers of pixels per subzone. Then, equation (TT~5T > rewrites 

(23) 3r-3^=\^lh 

2n h n h 1 

where nj is the number of cells in the selected zone Z n j, and nu and nj 2 are respectively the number of cells in the 
subzones Z^j and Z^j . Although the decrease of the objective function provided by a cutting is proportional to the 
square of the refinement indicator associated with this cutting, unfortunately the proportionality coefficient depends on 
the cutting itself. Nevertheless, for high values of a, the function x i — > x ("-x) * s aml0St f at m tne neighborhood of a/2, 
and the highest indicator almost yields the best decrease of the objective function, at least during the first iterations of the 
algorithm. Furthermore, there is no more interaction between the pixels during the minimization process. Hence, there 
is no need to compute the refinement indicators for all zones at each iteration, but rather, it is advisable to compute only 
the indicators related to the two most recent subzones, and to keep in memory the values for the unchanged zones. An 
optimized version of the refinement indicators algorithm dedicated to the diagonal models is described in |4|. 

For this first multidimensional application of the algorithm, we have simply defined the pseudo-sign of the fine vector 
gradient on each pixel by the sum of the signs for each component. Let sgn be the sign function on the field of real 
numbers (i.e., sgn(x) = — 1 when x < 0, sgn(x) = +1 when x > and sgn(0) = 0). Then, we define the sign of any triple 
by sgn(x,y,z) = sgn(sgn(x) +sgn(y) + sgn(z)). As developed in subsection !3.31 this is a very basic heuristics to deal with 
the multidimensional case, but it is sufficient for our present illustration purposes. 

In subsection [372] the choice for the refinement strategy defining the set of cuttings and the subset of those with the 
highest indicators was let to the user. Here, we consider two strategies. The first one is a mild flavor of the best strategy 
consisting in choosing at each iteration the cutting related to the highest indicator which guarantees in the case of the 
identity model the best descent of the objective function. In fact, we follow the pseudo-sign of the gradient as defined 
above. This strategy actually provides an image segmentation technique. The second one is called the dichotomy strategy. 
It is an example of strategy based on elementary cuttings deployed in all zones. The only cuttings allowed are those 
splitting a zone vertically or horizontally in its middle. This illustrates what could happen in the general case where the 
best strategy is not necessarily optimal and where the adaptive parameterization technique is based on the choice for a set 
of cuttings. Notice that this could provide interesting creative filters for image processing. 

The next Figures follow the same scheme. They are all made of four images. The data image d is displayed on the 
lower left corner. The segmented image y?° pt is displayed — in true colors — on the lower right corner, and the corresponding 
n-zone zonation T n is displayed in false colors on the upper right corner. The pseudo-sign of the fine gradient Vj? (pn Pt ) is 
displayed on the upper left corner where positive values are depicted in white, negative values in black and small absolute 
values in gray. Some monitoring counters are also displayed. This is actually the interactive display of the OCaml driver 
at runtime, see 0. 

The first series of experiments exploit the best strategy, hence we may expect the fastest decrease of the objective 
function. The results obtained after iterations number 1, 2, 5 and 10 are respectively represented in Figure [2][3]|4] and [5] 

At the beginning, in Figure|2l the initial zonation contains only one zone, the initial parameter is m 1 "" = (0,0,0) — i.e., 
a black image — and the initial value of the objective function is j[ mt = 8.0 10 9 . After the initialization phase, the objective 
function is j/ j 01 " = 1.2 10 9 , meaning that 62% of the data are explained^. Then, the pseudo-sign of the first fine gradient g° pt 
already depicts the head of the kangaroo. This cutting is associated with the — high — refinement indicator l^ est = 1.0 10 7 . 

After the fifth iteration, in Figure |U 90% of the data are explained and the main features of the original picture are 
visible. The last highest indicator was already only 7.4% of the first one. With 10 zones, in Figure[5] 94% of the data are 
explained. For the eye, the segmented image is very close to the original one. The last highest indicator was a flat 1 .7% 
of the first one. 



The second series of experiments use the dichotomy strategy. The result obtained after 2 and 100 iterations are 
respectively represented in Figures|6]and|7] 

The first iteration is exactly the same as with the best strategy in Figure[2] Then, the first cutting does not reveal at all 
the head of the kangaroo. The image with two zones, in Figure|6] only explain one percent more of the data than the one 



5 At iteration n, the percentage of data explained is defined as 100 x I 1 — 
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with one zone. The corresponding indicator is only 28% of the first indicator with the best strategy. Indeed, the chosen 
family of cuttings is not well suited for the problem. Even at iteration number 100, in Figure[7] only 87% of the data are 
explained and the last highest indicator was about 1 % of the first one. 

Of course, this strategy is inefficient for the segmentation of images, but exploring a subset of possible cuttings at each 
iteration acts as a preconditioner, and it may also help in avoiding to fall into local minima in the general case. Notice that 
after 100 iterations of the dichotomy strategy, the gradient keeps very sharp details that has already disappeared after only 
10'iterations of the best strategy. We could have done a better job by allowing any vertical and horizontal cuttings, but 
there would have been much more indicators to evaluate. The dichotomy run took 54 seconds, compared to the 3 seconds 
for the best run (only add tind and topt, as the total time ttot also takes into account the display). 

5 Conclusions 

We give a general formulation of the refinement indicator algorithm for the estimation of vector valued distributed param- 
eters in any partial differential equation. 
The main findings are: 

(i) in the multidimensional case, the refinement indicator associated with a new degree of freedom is the norm of the 
corresponding Lagrange multiplier. 

(ii) In the linear case, the refinement indicator also measures the decrease of the optimum least-squares data misfit objec- 
tive function when the corresponding degree of freedom is open. 

(iii) When applied to the inversion of the identity, the vector valued version of the refinement indicators algorithm provides 
a surprisingly powerful tool for the segmentation of color images (where colors are three-dimensional vectors in the RGB 
model). 
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Figure 2: Best strategy, 1 zone. Bottom: data image d (left), segmented image pj (right). Top: pseudo-sign of the 
gradient g° pt (left), 1-zone zonation f?i (right). 




Figure 3: Best strategy, 2 zones. Bottom: data image d (left), segmented image p 2 (right). Top: pseudo-sign of the 
gradient g 1 ^ 1 (left), 2-zone zonation 2>2 (right). The refinement indicator for the first cutting was /j best = 1 .0 10 7 . 
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iter = 5 nind = 10 nopt = 5 nzon = 5 

ttot = 11.9 tint! = 0.415 topt = 0.916 cost = 8.09e+07 




Figure 4: Best strategy, 5 zones. Bottom: data image d (left), segmented image p 5 (right). Top: pseudo-sign of the 
gradient g^ pt (left), 5-zone zonation 2>5 (right). The refinement indicator for the fourth cutting was l^ est — 7.5 10 5 . 




iter = 10 nind = 45 nopt = 10 nzon = 10 

ttot = 24.3 tind = 1.02 topt = 1.89 cost = 2.76e+07 




Figure 5: Best strategy, 10 zones. Bottom: data image d (left), segmented image (right). Top: pseudo-sign of the 
gradient g°^ 1 (left), 10-zone zonation fPio (right). The refinement indicator for the ninth cutting was 1^ cst = 1.7 10 5 . 
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iter = 2 nind = 2 nopt = 2 nzon = 2 

ttot = 4.07 tind = 0.152 topt = 0.392 cost = l.lle+09 




Figure 6: Dichotomy strategy, 2 zones. Bottom: data image d (left), segmented image p 2 p (right). Top: pseudo-sign of 
the gradient g ^ 1 (left). 2-zone zonation 2>2 (right)- The refinement indicator for the first cutting was / dlchotom y — 2.9 10 6 . 




Figure 7: Dichotomy strategy, 100 zones. Bottom: data image d (left), segmented image p,Q (right). Top: pseudo-sign 
of the gradient g ^ (left), 100-zone zonation !?ioo (right). The refinement indicator for the 99th cutting was J g9 c otomy — 
2.8 10 4 . 
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